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Steady-state thermodynamics of non-interacting transport be¬ 
yond weak coupling 
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Abstract - We investigate the thermodynamics of simple (non-interacting) transport models 
beyond the scope of weak coupling. For a single fermionic or bosonic level - tunnel-coupled to 
two reservoirs - exact expressions for the stationary matter and energy current are derived from 
the solutions of the Heisenberg equations of motion. The positivity of the steady-state entropy 
production rate is demonstrated explicitly. Finally, for a configuration in which particles are 
pumped upwards in chemical potential by a downward temperature gradient, we demonstrate 
that the thermodynamic efficiency of this process decreases when the coupling strength between 
system and reservoirs is increased, as a direct consequence of the loss of a tight coupling between 
energy and matter currents. 
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Introduction. — Thermodynamic studies on 
nanoscale processes have been investigated intensively 
iver the last years. For nanosystems, the traditional pic¬ 
ture of a Carnot engine that is connected in an alternating 
fashion to different thermal reservoirs is often replaced 
by stationary currents flowing through a quantum system 
continuously connecting such reservoirs [^. When the 
coupling between system and environments is sufficiently 
weak, a master equation approach is well applicable. 
Within standard techniques it is possible to derive 
thermodynamically consistent master equations of Lind- 
blad form, which lead to a positive entropy production 
and thus obey the laws of thermodynamics. 

With some exceptions [^, the thermodynamic study of 
such master equations is usually limited to the weak cou¬ 
pling regime between the system and reservoirs . Within 
such master equation descriptions, the energy levels of the 
system become infinitely sharp. In particular for a single 
intrinsic transition frequency, every particle that passes 
through the system must then carry a defined amount of 
energy, resulting in a direct proportionality of the matter 
and the energy flow (termed tight coupling). 

More recently, the fate of thermodynamic bounds 


has also been investigated beyond master equation ap¬ 
proaches [^[^. However, many studies were limited to the 
case of a single reservoir that does not support stationary 
currents [^. 

In this paper, we will consider stationary transport 
through either a bosonic or fermionic site between two 
reservoirs of similar type. The weak-coupling results are 
of course faithfully reproduced but for the particular prob¬ 
lems we will be mainly interested in the strong-coupling 
regime. 

The simplest transport model. — To perform our 
thermodynamic studies we investigate the simplest con¬ 
ceivable models of quantum transport. To support a sta¬ 
tionary current, such models should encompass (at least) 
two reservoirs, tunnel-coupled by a simple system that 
locally has only a single transition frequency. The Hamil¬ 
tonian is thus of the form 


H — ed^ d p ^ ^ 


kct 


+ ^ ^ > (1) 

koc 
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where e denotes the system-intrinsic transition frequency, 
ujka the frequency of mode k of reservoir a G {L,R}, and 
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tka the corresponding tunneling amplitudes between sys¬ 
tem and reservoir. The annihilation operators for the sys¬ 
tem (d) and the reservoirs {cka) can be of either bosonic 
(upper sign) or fermionic (lower sign) nature, manifested 
in corresponding commutation or anticommutation rela¬ 
tions, respectively. In the following expressions, we will 
adopt this convention, i.e., in case of differences between 
bosons and fermions the upper sign will hold for the 
bosonic transport model, and the lower for the fermionic 
counterpart. 

We will solve for the resulting global dynamics, using 
only that initially the reservoirs are prepared in grand- 
canonical equilibrium states, such that the initial state of 
the compound system is given by 

Po — Ps Tr i 

Qe{L.R} t / 

with Tr |/5g} = 1. Here, Pa and Ha denote the initial in¬ 
verse temperature and chemical potential of reservoir a, 
Ha its Hamiltonian - cf. Eq. § - and Na^Ek 4aCfca 
the corresponding particle number operator. Obviously, 
there is no initial entanglement between the system and 
the environment. Our treatment does not require the am¬ 
plitudes tka to be small. However, we remark that in 
general, for large amplitudes the distinction between sys¬ 
tem and reservoir becomes somewhat fuzzy, since some 
eigenstates of the Hamiltonian will extend over both sys¬ 
tem and reservoir. For example, it is in principle possi¬ 
ble that energy contained in the interaction Hamiltonians 
then contributes non-negligibly to steady-state transport. 
We will see later that for our setup this is not the case, 
enforced by the initial condition and the assumption 
of a continuous spectral coupling density. 

Equation-of-motion method. — In the equation-of- 
motion-method, one simply solves the Heisenberg picture 
dynamics = i [H, A] for the bosonic or fermionic op¬ 
erators generated by Eq. Q exactly. Afterwards, physi¬ 
cal observables of interest can be obtained from the time- 
evolved operators. 

For a quadratic Hamiltonian without interactions as 
considered here, the Heisenberg equations of motion just 
yield a linear set of first order differential equations 

d = ied T i ^ ^ tka^ka J 

ka 

^ka — ^^ka^ka ^ ^^ka^ 7 (3) 

and similarly for the creation operators, where the bosonic 
(—) or fermionic (-I-) nature of the operators just manifests 
in the sign. We note that in case of interactions present, 
the set will not close and approximations must be applied 
for further treatment. Similar to coupled equations for the 
evolution of amplitudes [^, this closed set of infinitely 
many differential equations can be transformed to an al¬ 
gebraic set of equations using a Laplace transform 


C{d} = d{z) = /y e~^*d{t)dt, which can then be solved 
explicitly in terms of the initial operators (denoted by d 
and Cka, respectively). For later convenience, we state the 
result for d{z) and 4a (^) explicitly 


d{z) = 


^ + + 
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( 4 ) 


Now, exact solutions can in principle be obtained by in¬ 
verting the Laplace transform, which in the general case is 
typically tedious. Since we want to consider steady state 
configurations, we must assume an infinite and continu¬ 
ous distribution of reservoir modes k (finite-sized quantum 
systems evolve periodically). Technically, this corresponds 
to the replacement of 
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\tkX 
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with the energy-dependent tunneling rate (or spectral cou¬ 
pling density) rc(w) = 2'KEk Eal S{uj - ujka)- Finally, 
we mention that to obtain the solution for the hermi- 
tian conjugate operators, the Laplace transform variable 
should be kept invariant to remain consistent with our fur¬ 
ther calculations (meaning the same Laplace transform is 
applied to creation and annihilation operators in the time 
domain). 


The currents. — We are interested in the long-term 
matter and energy currents entering the right lead and 
therefore compute the Heisenberg equation of motions for 
the particle number operator and the Hamiltonian of the 
right reservoir, yielding 


-^R — T1 ^ ^ |^^fcR4R,^ h.c. , 

k 

— 4“1 h.c. 

k 


( 6 ) 


Computing expectation values and considering the long¬ 
time limit, we find that the steady-state matter and en- 
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ergy currents into the right reservoir can be calculated by 
convolution integrals in Laplace space 
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lim (-/Vr) = lim zL 

t^oo 2—>-0 
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CT — ie + C'L(d) + C'R(cr) = 0 can only be solved by a purely 
imaginary a = iay with (for e > 0) CTy > 0. Since we have 
that lim [z — cr + ie + C'l(z — a) + Cr(z — cr)]^^- = 0, 
invoke the rule of L’Hospital to 


we 

find 


“ cr))d(T + h.c., 


can 
that 
1 


z — a+i€+Ci,{z — a)+Cn{z—<j) 

= 0 , where we have used that 


(Ty > 0 leads to a divergence of the integral. Therefore, 
as expected, the steady-state currents will not depend on 
the initial occupation of the system. 

The same arguments can be applied to evaluate the 
other contributions. In particular, combining the first 
term in c^.r(o') with the second term in d{z — a) we obtain 


( 7 ) a stationary contribution to the matter current 


Here, the real constant 7 has to be chosen such that the 
poles of c|,j^(tT) lie left to it and the remaining poles of 
d{z — a) are on the right of the integration contour as 
depicted in Fig. 
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for which we used 



Figure 1: (Color Online) Sketch of the relevant poles arising in 
the expectation value {cl-^(a)d{z — a)). The contour parameter 
7 must be chosen such that the poles of clj^(a) (blue) lie left 
and the poles of d(z — a) (red) lie right to it. When the limit 
2 —> 0 is performed afterwards (thin dashed arrows), it becomes 
visible that in order to yield a long-term contribution, poles 
inside the integration contour described by cr — ica = 0 can 
only contribute if they have a conjugate pole (z — a + iu = 0) 
outside the contour. In this case, when evaluating the residue 
one obtains lim —= 1 . 

- vn 2 —CT + ICl) 

Z —I c 7 = 1 Ll) 

To evaluate the currents explicitly, we note that 
due to our condition of initial thermal equilibrium 
only few terms may potentially contribute, since 
((^k'a'Cka) = ^fcfe'dcta'Ra(wfe) with the Bose (—) and Fermi 
(-I-) distributions 


lim Caiz — ica) = C'o,(—ica) 
2—>-0+ 


2 


+ iSQ(ca), 


( 10 ) 


where ^^(ca) = dca' (with Cauchy principal 

value V) represents a coupling-induced level renormaliza¬ 
tion (often called Lamb shift). Combining the last term 
in c])p^((t) with the second term in d{z — a) one obtains a 
second contribution 


r(2) 


00 


0 


doj CR(ii:a) [rL(a;)nL(ca) + rR(u;)nR(ca)] 


( 11 ) 


where F(ca) = ^^(w) and in a similar fashion for S(ca). 

To further separate real and imaginary parts we can use 
lim C'R(cr) = C'R(ita) = - iSR(ca). 

(T—vicj 

Combining these expressions in the total matter current 
Im = 23fJ(/M^ + and performing similar calculations 
for the energy current, we finally obtain for the steady- 
state currents the expressions 


Bm 


Ie 


00 

J Gc{u}) [nL(w) - nR(a;)] 5'c(a;)da;, 

0 


to ■ Gc(w) [nL(w) - nR(w)] S'c(a;)da;, (12) 


-Me) =p 1 ■ 

We can evaluate these contributions with the residue the¬ 
orem. 

For simplicity, we first discuss the term arising from 
the combination (d^d) = uq. Here, the residue arising 
from the pole at cr = -l-iaifco, is multiplied by z and 
thus vanishes as z —>■ 0. Furthermore, the equation 


with the factors 


Gc(w) 


5c (w) 


rL(w)FR(a;) 

’ 

1 r(,.)/2 

’<{r{u)i2Y + \u-r.-n<j)f’ 


and where r(a;) = rL(w) -I- rR(tLi) and E(w) = El(w) -I- 
I]R(a;). The difference between bosonic and fermionic 
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transport enters in the different distributions Q and the 
fact that for bosons, the chemical potentials must be neg¬ 
ative (otherwise the total particle number does not con¬ 
verge), whereas they are unbounded for fermions. Further¬ 
more, also for fermions we assumed the reservoir frequen¬ 
cies uJka to be positive. Relaxing this assumption would 
simply extend the lower bound in all integrals to — oo. 

The Landauer form of the current is well-known for 


purely fermionic 12 -15 and purely bosonic transport 


16 


Structurally similar results hold for bosonic transport 


through a central spin 18 19 


A nice feature of Eqns. (12) and (13) is the common rep¬ 


resentation for both bosons and fermions. The key advan¬ 
tage of this Landauer-\ike representation 20-22 is that it 


iM 
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rL(e)+rR(e) 
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(15) 


rents (12) can be interpreted as a broadening of the system 


energy level, since all modes contribute to the transport, 
weighted by the tunneling rates. This broadening violates 
the tight-coupling of matter and energy fluxes. 

We remark that the long-term energy content of the in¬ 
teraction Hamiltonian between system and the right lead 
can be extracted fro m t he imaginary parts of the integrals 
in Eq. ^ and Eq. (U) via |(iTi,R(t))| -t 23(/^^ 

Since these remain finite for large times, we conclude that 
the interaction energy does not contribute to the steady- 
state energy currents. For finite-time statements how- 
these contributions will matter. 


ever 


23 26 


Thermodynamics. — 

Steady-state entropy production. When a system is 
coupled to different equilibrium environments, its change 
of entropy can be split in terms that have an interpretation 


as entropy flow and entropy production [^. In particular, 
when a hnite system reaches a steady-state (and hence its 
change of entropy vanishes), the internal entropy produc¬ 
tion rate A^S must be balanced by the entropy flow Ae.S' 
entering the system. 


AiS = -A,S = -J2 PcQo. > 0. 


(16) 


Here, Qa denotes the heat current entering the system 
from reservoir a. We note that positivity of integrated en¬ 
tropy production - defined as difference between change 
of internal entropy and integrated entropy flow - has been 
proven generally and also at the level of individual 


is completely independent of the actual shape of rQ(a;). 

For consistency, we note that in the weak-coupling 
limit rQ,(w) —)■ 0 (also implying T.{uj) —>■ 0), one of 
the factors becomes a Dirac-Delta function S'c(w) —t 
5{oj — e — E(a;)) —t (5(a; — e), such that the integrals col¬ 
lapse and the master equation results are reproduced 


trajectories 28 . Here, we will demonstrate explicitly that 

hout changes in the basic ther- 
- at steady state supports a 


our transport scenario - without changes in the basic ther¬ 
modynamic definitions 
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Assuming both strong coupling and flat tunneling rates 
ra(a;) —>■ Tq [such that S'c(w) —t 2/(7rr)], we can obtain 
analytic results for the bosonic (upper sign) and fermionic 
(lower sign) matter and energy currents 


positive entropy production rate (161. The proof is simple 

Using 


and general for Landauer-Biittiker transport [29| . 
matter and energy conservation (one can of course explic¬ 
itly check that at steady state the first law is obeyed), the 
heat currents entering the system are defined as 

Qe = +4e — IJ-eIm , Qr = —Ie + P-kIm ■ (17) 

Altogether, the internal entropy production rate is thus of 
the form 


= (/3r — Pe)Ie + {pePe — PkPk)Im ■ 


(18) 


Inserting the integral expressions for the currents to the 
equation above yields 


AiS' = / x(a;)Gc(w)S'c(w)da;, 


(19) 


where a = ^ (rL+ra)^ Li 2 (a;) denotes the polylog func¬ 
tion. Even simpler expressions arise for fermions when the 
integral is extended to the complete real axis (not shown). 

For finite coupling strengths, the factor Sc{aj) encodes 
the coupling-induced level renormalization via E(a;). The 
integral over all energies w in the expressions for the cur- 


X(w) = [w(/3r - /3l) + (pePe - PrPr)] [nL(w) - nR(a;)] . 

As one obviously has Gc(w)<S'c(w) > 0 Vw, since 
r(w) > 0, it is sufficient to show that x(a;) > 0. In the 
trivial case of equal temperatures and chemical potentials 
(nL(w) = nR(a;)), all currents vanish and the entropy pro¬ 
duction is zero. Only for different temperatures, a non¬ 
trivial root of xi^) <^0 = exists with 




= 0 , 


do; 2 


X(w) 


(Pe - PkY 


cosh 




( 20 ) 

> 0 . 


Tl 


Consequently, the root at wq also corresponds to a global 
minimum. As x^(<^) is continuously differentiable for fi¬ 
nite temperatures and vanishes in the limit of a; —>■ oo, we 
obtain that x(w) > 0. Non-negativity of the integrand in 
Eq. (191 along the real axis thus also implies a positive 


internal entropy production rate at steady state. A sim¬ 
ilar proof would hold if the Bose-Einstein or Fermi-Dirac 
distributions were replaced by Boltzmann factors. 
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We mention that this second law automatically implies 
that at equal temperatures matter will always flow from 
reservoirs with large to the one with small chemical poten¬ 
tial and that at equal chemical potentials heat will always 
flow from hot to cold reservoirs. 


Nanothermal engine. Now we consider a configuration 
in which our transport setup may function as a simple 
nanothermal engine. For fermions, exactly this kind of 
model has been considered before in the weak-coupling 
Assuming that the temperature of the left 


30 


regime 

lead is smaller than that of the right /3 l > /3r, whereas the 
potential gradient V = — /j,r > 0 is tilted in the other 

direction, it is in certain parameter regions possible to 
induce a current against the potential gradient (Im < 0), 
such that the generated power 


F = —/m(ml —Mr) (21) 


becomes positive. For positive power, the thermodynamic 
efficiency of this process is defined as the ratio between 
generated power and the heat current entering the system 
from the hot (right) reservoir, 


V = 


1^m(ml — Mr) 
Ie — MR-^m 


0(-1^m(ml - Mr)) : (22) 


where we have introduced the Heaviside-function to en¬ 
force the correct regime (finite-time realizations [26|[M][^ 


would have to relate work output and heat input). From 
positivity of Eq. ( [I^ we can deduce that the efficiency 
is always upper-bounded by Carnot efficiency g < rjc = 
1 - Tl/Tr. 

Loss of tight-coupling. In the weak-coupling limit - 


compare Eq. (14) - there is only a single transition fre¬ 


quency e in the system, and the currents become tightly 
coupled /e = c/m. Then, one can explicitly show that 
Carnot efficiency is actually reached when the current (and 
hence, also the power) vanishes. Due to the fact that 
the power at maximum efficiency is zero, it has become 
customary to consider the efhciency at maximum power 
output instead. Since this still goes beyond the linear re¬ 
sponse regime, this becomes an optimization problem that 
can in general not be solved analytically but requires nu¬ 
merical approaches. 

In contrast, in the strong-coupling regime - compare 
Eq. (15) ~ one can numerically verify that the power 


at maximum efficiency does not need to be significantly 
smaller than the maximum power. 

As we will show in the following, beyond the weak- 
coupling regime, the power generated at maximum effi¬ 
ciency remains finite. With the aim of a symmetric and 
consistent description for both bosons and fermions, we 
introduce the new variables 


Thus, with the correct choice of the parameters /3 and 
M, we are able to vary the (inverse) temperature gradient 
A/3 and the voltage bias Am over the complete real axis. 
In particular for the bosonic model, this representation is 
useful to ensurse negative chemical potentials throughout. 
In Fig. we plot the efficiency versus the temperature and 
potential gradients in the weak (left) and strong (middle) 
coupling regimes and also provide a direct comparison of 
the different coupling strengths (right). We note that in 
contrast to the master equation efficiency - cf. Eqns. (221 
and (14) - the true efficiency also varies with respect to 
temperature. It is visible that the maximum efficiency 
versus bias voltage will decrease with increasing coupling 
strength. Whereas for small coupling (left) it is close to the 
master equation results and thus close to Carnot efficiency, 
it decreases drastically for larger couplings (middle). We 
attribute this to the loss of tight-coupling in the strong¬ 
coupling regime. As a more positive feature of the strong¬ 
coupling regime - in particular for practical applications 
- however, we also observe that the curves of efficiency 
at maximum power and maximum efficiency are closer to¬ 
gether and that maximizing the efficiency does not require 
to accept zero power. Unfortunately, also in the strong¬ 
coupling regime - compare Eqns. (151 - these curves will 
not coincide, such that one still has find a tradeoff between 
maximum power and maximum efficiency 33 . The plots 


for bosonic transport (not shown) are very similar. 


Summary. — Using a simple equation of motion tech¬ 
nique in combination with Laplace transform methods, we 
have calculated the exact steady-state energy and mat¬ 
ter currents of simple two-terminal transport models both 
for fermions and bosons, yielding a Landauer-type rep¬ 
resentation. The approach can be expected significantly 
more complicated both from a technical and conceptual 
perspective when applied for finite times. When com¬ 
bining the resulting expressions into heat currents, we 
showed explicitly that the steady-state entropy produc¬ 
tion rate is positive. Using this second law inequality, we 
discussed the efficiency of converting a thermal gradient 
into power. For our transport models we noticed that 
for increased system-reservoir coupling strength, the max¬ 
imum efficiency decreases. Importantly, we also note that 
the power at maximum efficiency is finite, rendering the 
strong-coupling regime a potential candidate for practical 
applications. 
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Figure 2: (Color Online) Plot of the renormalized fermionic heat-to-power-conversion efficiency rj/rjc versus dimensionless bias 
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efficiency and thus significantly smaller than Carnot efficiency. For strong coupling (/fF = 10.00, middle panel) efficiency at 
maximum power (dashed red) and maximum efficiency (dashed orange) are much closer, and the master equation curve for zero 
power (solid orange) now yields a wrong bound for the region of finite efficiency. The right panel shows for fixed temperature 
difference A/1//3 = 1.7 (top, horizontal dashed line in contour plots) and for fixed potential difference PAfj, — 1.75 (bottom, 
vertical dashed line in contour plots) that for increasing coupling strength the maximum efficiency decreases. Furthermore, 
it can be seen that also for small coupling strengths, the exact solution predicts that the maximum efficiency is reached 
slightly before the power vanishes (magenta). Contours denote efficiency steps of 0.05. Symmetric Lorentzian tunneling rates 
rL(tu) = rR(tj) = + S^) with S = e have been used throughout. Other parameters: /3e = +2 and = —e/2. 
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